Phylo Tree trSusD Prevotella

Last modification: 20 May, 2020

@EricJCGalvez

# check.packages function: install and load multiple R packages.
# Check to see if packages are installed. Install them if they are not, then load them into the R session.
check.packages <- function(pkg){
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    if (length(new.pkg)) 
        install.packages(new.pkg, dependencies = TRUE)
    sapply(pkg, require, character.only = TRUE)
}


# Usage example
packages<-c("tidyverse", "ggplot2", "ggthemes", "knitr", "DT", "ggpubr", "ggsci", 
            "data.table", "ggtree", "ape", "treeio", "here")
check.packages(packages)

load tree

tree <- treeio::read.newick("../data/RAxML_bipartitions.susD_rooted.newick")

##Checking the clade number
## Is necesary to plor first removing the " aes(color=group)"
## then running the command MRCA and the plot interval is posible to get the node number

#node_rnode_ref <-252
tree2 <- groupClade(tree, .node = c(181, 221, 230))


df <- as.tbl(read_delim("../data/tree_io_metadata_susD_R.txt", delim = "\t"))

## Parsed with column specification:
## cols(
##   tip_id = col_character(),
##   Genome_MAG = col_character(),
##   TipLabels = col_character()
## )

## Warning: 1 parsing failure.
## row col  expected    actual                                  file
## 150  -- 3 columns 2 columns '../data/tree_io_metadata_susD_R.txt'

rownames(df) <- df$tip_id               

## Warning: Setting row names on a tibble is deprecated.

setDT(df, keep.rownames = TRUE)[]

##                           rn                  tip_id       Genome_MAG
##   1: BACOVA_03427_XylL_SusD1 BACOVA_03427_XylL_SusD1       BACOVA_Xul
##   2: BACOVA_03429_XylL_SusD2 BACOVA_03429_XylL_SusD2       BACOVA_Xul
##   3:  BACOVA_04392_XylS_SusD  BACOVA_04392_XylS_SusD       BACOVA_Xul
##   4:            MAG609_00227            MAG609_00227           MAG609
##   5:            MAG609_00350            MAG609_00350           MAG609
##  ---                                                                 
## 146:              PROD_01883              PROD_01883 Mouse Prevotella
## 147:              PROD_02119              PROD_02119 Mouse Prevotella
## 148:              PROD_02217              PROD_02217 Mouse Prevotella
## 149:              PROD_02444              PROD_02444 Mouse Prevotella
## 150:              PROD_02831              PROD_02831 Mouse Prevotella
##         TipLabels
##   1: BACOVA_03427
##   2: BACOVA_03429
##   3: BACOVA_04392
##   4:         <NA>
##   5:         <NA>
##  ---             
## 146:         <NA>
## 147:         <NA>
## 148:         <NA>
## 149:         <NA>
## 150:         <NA>

colnames(df)[1] <- "label"

tree_df <-  treeio::full_join(tree2, df, by = 'label')


###### To check node number ID ## 
# susc1_xylL !

#susd2_xylL <- MRCA(p, "BACOVA_03427_XylL_SusD1", "MAG612_00689") #230
#susd1_xylL <- MRCA(p, "BACOVA_03429_XylL_SusD2", "MAG609_01227") #221
#susd1_xylS <- MRCA(p, "MAG611_01469", "MAG612_00422") #181


tree_df@data$group

##   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0
##  [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 0 0
##  [75] 0 0 0 0 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0
## [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2
## [223] 2 0 0 0 0 0 0 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [297] 0 0 0
## Levels: 0 1 2 3

levels(factor(df$Genome_MAG))

## [1] "BACOVA_Xul"       "MAG609"           "Mouse Prevotella" "P. copri MAG"    
## [5] "PINT"

my_col_host= c("0"="black", "1"="#2166ac", "2"="#2166ac", "3"="#2166ac",
               "BACOVA_Xul"="#fec44f","MAG609"="#cb181d",
               "Mouse Prevotella"="#a8ddb5", "PINT"="#fb8072", "P. copri MAG"="#525252")

p <- p <- ggtree(tr = tree_df, layout = "rectangular", aes(color=group), right = TRUE , size=1) +
      geom_tippoint(aes(color=Genome_MAG), size=1.5, alpha= 0.7) +
      scale_color_manual(values=my_col_host) +
      geom_rootedge(0.01)
p

## number infered from MRCA function
clade1 <- 181
clade2 <- 221
clade3 <- 230

p1 <- scaleClade(p, node=252, scale=1) %>% 
  scaleClade(., node=243, scale=1)  %>% 
  scaleClade(., node=225, scale=1)

p2 <- p + geom_hilight(node=clade1, fill="#feb24c", alpha=.3, extend = 1) +
    geom_hilight(node=clade2, fill="#feb24c", alpha=.3, extend = 1) +
    geom_hilight(node=clade3, fill="#feb24c", alpha=.3, extend = 1) +
    geom_tiplab(aes(label=TipLabels, na.omit=T), linetype=NA, align = F,
                  vjust=0.5, hjust=-.1, color= "black",
                  size=1.5, offset = .005)

p2

#export 6x9

p <- ggtree(tr = tree_df, layout = "fan", aes(color=group, fill=Genome_MAG), size=1,  #0.3
            open.angle = 15, right = T) +
      #geom_tippoint(aes(color=Genome_MAG), size=1, alpha= 0.7, shape= c(10)) +
      geom_point(aes(shape=Genome_MAG, color=Genome_MAG, alpha= 0.9, size=2), na.rm=TRUE) +
      scale_shape_manual(values = c(16, 20, 20, 20, 1)) +
      scale_color_manual(values=my_col_host) +
      scale_fill_manual(values=my_col_host) +
      geom_rootedge(0.5)
p

p1 <- scaleClade(p, node=181, scale=1) %>% 
      scaleClade(., node=221, scale=1)  %>% 
      scaleClade(., node=230, scale=1)


clade1 <- 181
clade2 <- 221
clade3 <- 230

p3 <- p1 +
    geom_hilight(node=clade1, fill="#feb24c", alpha=.3, extend = 1) +
    geom_hilight(node=clade2, fill="#feb24c", alpha=.3, extend = 1) +
    geom_hilight(node=clade3, fill="#feb24c", alpha=.3, extend = 1) +
    geom_tiplab(aes(label=TipLabels, na.omit=T), linetype=NA, align = F,
                  vjust=0.5, hjust=-.1, color= "black",
                  size= 2, offset = .005) 

p3

rotate_tree(p3, 180) 

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#export 6 x 9